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1  Summary  of  Principal  Accomplishments 

The  research  conducted  over  the  past  four  years  of  this  AFOSR  grant  has  focused  on  one 
of  the  fundamental  factors  influencing  the  accuracy  of  floating-point  computations,  namely, 
the  condition  or  sensitivity  of  the  particular  problem  being  solved.  Many  fundamental 
calculations  in  scientific  and  engineering  computation  in  general,  and  control  engineering  in 
particular,  can  be  viewed  as  functions  that  map  a  set  of  input  values  to  a  set  of  output  values. 
Assuming  the  algorithm  used  in  the  computation  is  numerically  stable,  the  sensitivity  of 
this  input-output  map  determines  how  accurately  the  output  values  can  be  estimated.  The 
methods  used  in  our  research  rely  on  the  fundamental  results  in  [1]  and  [19]  wherein  it  is 
shown  that  significant  information  about  large  linear  operators  (which  arise  naturally  as 
the  derivative  of  some  nonlinear  function  of  interest)  can  be  obtained  from  their  action  on 
smaller-dimensional  subspaces.  Such  information  can  be  determined  reliably  and  efficiently 
using  the  small-sample  statistical  method.  In  addition  to  further  advancing  our  basic  theory, 
we  have  explored  applications  of  statistical  techniques  to  a  wide  variety  of  problems  in 
control  theory  and  image  processing.  Our  research  has  been  reported  in  over  40  scholarly 
articles. 

1.1  Sensitivity  in  Linear  Algebra  and  Control  Computations 

Standard  approaches  to  measuring  the  condition  of  various  problems,  such  as  solving  a. 
system  of  linear  equations,  compress  all  sensitivity  information  into  a  single  "condition 
number."  Thus,  a  loss  of  information  can  occur  in  situations  in  which  this  standard  con¬ 
dition  number  does  not  accurately  reflect  the  actual  sensitivity  of  a  solution  or  particular 
entries  of  a  solution.  Our  new  method  overcomes  these  and  other  common  deficiencies  by 
measuring  the  effects  on  the  solution  of  small  random  changes  in  the  input  data.  By  properly 
scaling  the  results,  we  can  obtain  condition  estimates  for  each  entry  of  a  computed  solution. 
This  approach,  which  is  referred  to  as  small-sample  statistical  condition  estimation  (SCE), 
applies  to  both  linear  and  nonlinear  problems.  In  the  former  case,  when  an  explicit  Frechet 
derivative  of  the  computed  quantity  in  question  is  available,  the  method  is  especially  effi¬ 
cient,  costing  no  more  than  standard  normwise  or  componentwise  estimates.  Moreover,  SCE 
has  the  advantage  of  considerable  flexibility.  For  example,  it  easily  accommodates  restric¬ 
tions  on.  or  structure  associated  with,  allowable  perturbations.  The  method  has  a  rigorous 
statistical  theory  available  for  the  probability  of  accuracy  of  the  condition  estimates.  We 
have  applied  SCE  to  several  basic  problems  in  numerical  linear  algebra  and  several  papers 
have  been  published  or  are  in  press  in  the  leading  journals,  including  papers  on  solving  lin¬ 
ear  equations  [40],  linear  least-  squares  problems  [41].  and  eigenvalue/eigenvector  problems 
[37].  We  have  also  made  considerable  progress  in  applying  SCE  to  various  computational 
problems  in  control,  especially  the  numerical  solution  of  algebraic  Riccati  equations  [11], 
[15].  [30].  [31]. 

A  recent  thrust  in  our  research  program  has  been  devoted  to  new  algorithms  for  eval¬ 
uating  certain  key  matrix  functions  of  interest  to  control  engineers.  Perhaps  surprisingly 
to  some,  much  remains  to  be  done  before  algorithms  for  such  matrix  functions  as  the  ex¬ 
ponential  and  logarithm,  as  implemented  in  state-of-the-art  codes  in  Matlab’s  expm  and 
logm  for  example,  can  be  said  to  be  truly  reliable.  We  have  explored  a  so-called  Schur- 
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Frechet  method  [42]  of  evaluating  matrix  functions  which  consists  of  putting  the  matrix 
in  upper  triangular  form,  computing  the  scalar  function  values  along  the  main  diagonal, 
and  then  using  the  Frechet  derivative  of  the  function  to  evaluate  the  upper  diagonals.  This 
approach  requires  a  reliable  method  of  computing  the  Frechet  derivative  accurately.  For 
the  logarithm  this  can  be  done  by  using  repeated  square  roots  and  a  hyperbolic  tangent 
form  of  the  logarithmic  Frechet  derivative.  Pade  approximations  of  the  hyperbolic  tangent 
lead  to  a  Schur- Frechet  algorithm  for  the  logarithm  that  avoids  problems  associated  with 
the  standard  "inverse  scaling  and  squaring"  method.  Inverting  the  order  of  evaluation  in 
the  logarithmic  Frechet  derivative  gives  a  method  of  evaluating  the  derivative  of  the  ex¬ 
ponential.  The  resulting  Schur-Frechet  algorithm  for  the  exponential  gives  superior  results 
compared  to  standard  methods  on  a  set  of  test  problems  from  the  literature. 

Our  SC'E  method  can  be  applied  to  estimate  the  sensitivity  of  each  entry  in  a  matrix 
function.  For  example,  suppose  that  we  are  interested  in  the  sensitivity  of  the  (i.j)  entry 
of  a  matrix  function  F  =  F{A).  This  sensitivity  can  be  measured  by  the  norm  of  the 
gradient  of  the  map  .4  i —  Fij(A).  Using  SCE,  we  determine  sensitivity  by  perturbing  the 
argument  (,4)  in  a  random  fashion  and  looking  at  the  resulting  effect,  on  the  function  value. 
In  general,  a  random  perturbation  is  not  going  to  point  in  the  direction  of  the  gradient 
(i.e..  in  the  matrix  direction  that  produces  the  maximal  perturbation  in  Fiji  A)).  Because 
of  this,  a  scaling  factor  also  called  the  "Wallis"  factor,  must  be  introduced.  The  Wallis 
factor  depends  only  on  the  number  in  of  arguments  being  perturbed.  In  the  case  considered 
here,  m  =  n2  where  n  is  the  order  of  the  matrix  .4.  For  computational  purposes  the 
approximation  ^ ~  v/‘2/(  tt ( m  —  0.5))  is  sufficiently  accurate. 

This  condition  estimation  procedure  has  a  particularly  simple  form  when  expressed  in 
terms  of  a.  known  Frechet  derivative  and,  in  the  case  of  the  matrix  logarithm  and  matrix 
exponential,  such  a  derivative  is  known  explicitly.  This  random  perturbation  method  of 
condition  estimation  is  very  flexible  and  can  easily  be  adapted  to  problems  in  which  the 
perturbation  must  preserve  some  structure  of  the  matrix  such  as  symmetry. 

1.2  Sensitivity  in  Image  Processing 

During  the  past  few  years  we  have  worked  with  Dr.  Gary  Hewer  at  NAWC  in  China  Lake. 
CA  and  Prof.  B.S.  Manjunath  (UCSB)  on  sensitivity  issues  in  image  processing,  especially 
sensitivity  of  computed  results  for  segmentation  and  recognition  of  small  targets  in  im¬ 
ages.  registration  of  objects  in  sequences  of  images,  and  estimates  of  optical  flow  velocity 
fields.  The  goal  has  been  to  devise  procedures  such  as  the  SCE  method  for  automatically 
detecting  sensitivity  combined  with  efforts  to  develop  algorithms  for  feature  extraction  that 
are  inherently  less  sensitive.  A  general  framework  for  variational  approaches  to  image  seg¬ 
mentation  has  been  developed  that  includes  most  of  the  commonly  used  methods  such  as 
Mumford-Shah  and  Genian  type  functionals.  The  main  advantage  of  our  new  framework 
is  that  the  boundary  function  associated  with  the  segmentation  is  given  explicitly  in  terms 
of  the  approximation  function.  This  reduces  the  computational  burden  of  finding  the  seg¬ 
mentation  as  well  as  allowing  direct  analytical  comparisons  with  other  methods  with  regard 
to  sensitivity  to  noise  in  the  input.  A  description  of  the  variational  boundary  approach  is 
given  in  [32],  [43]. 

We  have  also  developed  a  second  approach  to  segmentation  and  edge  detection  based 


on  the  idea  of  "peer  group"  processing.  In  this  approach,  at  each  pixel  in  the  image  a  peer 
group  is  identified  based  on  nearness  of  intensity.  The  peer  group  can  also  be  based  on 
other  measures  such  as  nearness  in  texture  space  or  (for  optical  flow)  nearness  in  velocity 
space.  This  method  avoids  the  edge-blurring  difficulties  associated  with  gradient-based 
partial  differential  equation  approaches  that  do  not  distinguish  between  pixel  information 
for  different  regions  near  image  boundaries.  An  interesting  connection  between  the  shock 
front  image  processing  method  of  Rudin  and  Osher  and  the  peer  group  image  processing 
method  was  pointed  out  by  Stan  Osher  at  the  September  1996  OXR  workshop  at  UCLA. 
This  connection  allows  one  to  analyze  the  sensitivity  of  the  peer  group  method  with  respect 
to  additive  noise  and  to  conclude  that  the  total  variation  of  the  original  image  is  preserved. 
This  work  is  reported  in  [29].  [33].  [39]. 

1.3  Related  Topics 

The  Publications  list  below  includes  numerous  other  AFOSR-supported  papers  on  closely 
related  research  topics.  These  are  listed  below  and  were  discussed  in  detail  in  previous 
years*  Annual  Progress  Reports. 
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4  Interactioiis/Transitions 

4.1  Participation/presentations  at  meetings,  conferences,  seminars,  etc. 

During  the  tenure  of  this  grant,  the  P.I.  has  been  invited  to  deliver  three  major  plenary 
addresses  at  prestigious  control  theory,  numerical  linear  algebra,  and  applied  mathematics 
conferences.  The  first  was  at  the  Allerton  Conference:  the  second  was  at  the  Householder 
Symposium,  the  premier  venue  for  numerical  linear  algebra  and  held  only  every  three  years: 
and  the  third  was  at  the  recent  SIAM  National  Meeting,  held  last  year  at  Stanford  Lniver- 
sity.  Details  follow. 

T.  The  Matrix  Sign  Function  and  Its  Applications.  Invited  Plenary  Lecture  for  the 
Thirty-Third  Annual  Allerton  Conference  on  Communications,  Control,  and  Com¬ 
puting.  Allerton.  Illinois.  Oct.  4-6.  1995. 

2.  Statistical  Condition  Estimation  in  Numerical  Linear  Algebra.  Invited  Plenary  Lec¬ 
ture  for  the  Thirteenth  Householder  Symposium  on  Numerical  Linear  Algebra.  Pon- 
tresina.  Switzerland.  June  17-21.  1996. 

3.  Algorithms  for  Computing  Matrix  Logarithms  and  Exponentials.  Invited  Plenary 
Lecture  for  SIAM'S  45th  Anniversary  Meeting,  Stanford.  California,  July  14— lb. 
1997. 

4.2  Transitions  and  other  consultative  and  advisory  functions 

Prof.  Laub  consults  regularly  with  colleagues  at  The  MathWorks,  Inc.  in  Natick.  Mass, 
concerning  various  aspects  of  next-generation  software  and  algorithms  for  control.  Piof. 
Laub's  expertise  in  the  numerical  solution  of  algebraic  Riccati  equations  was  also  sought 
in  May  and  June  of  1996  by  Mr.  Kevin  J.  Shortelle  (sdi@afn.org)  of  System  Dynamics 
International.  Inc.  (9120  SW  46th  Blvd..  Gainesville,  FL  32608).  Shortelle  was  performing 
work  under  a  six-month.  880 Iv  contract  with  the  U.S.  Air  Force  to  investigate  the  bene¬ 
fits  of  innovative  nonlinear  estimation  techniques  for  improving  the  accuracy  of  integrated 
navigation  systems  (e.g.,  integrated  GPS.  Doppler,  and  INS  systems).  He  was  assessing 
not  only  conventional  nonlinear  approaches  (extended  Kalman  filter  (EKF),  iterated  EKF, 
modified-gain  EKF)  but  also  a  new  technique  referred  to  as  the  State  Dependent  Riccati 
Equation  (SDRE).  The  SDRE  techique  is  based  on  work  originating  from  research  at  the 
U.S.  Air  Force's  Wright  Laboratory,  at  Eglin  AFB.  Florida.  Dr.  Jim  Cloutier  (a  Wright 
Lab  scientist)  and  LTC  Curtis  Mracek  have  authored  several  papers  on  the  subject.  One  of 
Shortelle's  software  tasks  is  the  efficient  computation  of  solutions  to  the  ARE  based  on  cer¬ 
tain  system  state  matrices.  He  sought  from  Prof.  Laub  a  public-domain  Fortran  source  code 
available  for  such  a  computation  and  it  was  provided  to  him  together  with  some  explanation 
of  how  to  use  it  and  where  to  acquire  related  codes. 

Dr.  Kenney  has  had  conversations  and  provided  technical  guidance  concerning  condi¬ 
tioning  and  statistical  condition  estimation  for  automatic  target  recognition  and  image  pro¬ 
cessing  (such  as  optical  flow  computations)  with  Jim  Huang  (jim.huang@alphatech.com) 
of  Alphatech,  Inc.,  (50  Ball  Road.  Burlington.  MA  01803).  In  addition  to  discussing  sensi¬ 
tivity  of  optical  flow  algorithms,  Dr.  Kenney  sent  him  a  copy  of  a  Matlab  preprocessing 
routine  for  optical  flow  that  eliminates  a  lot  of  image  noise  from  saturated  or  dead  pixels. 


8 


5  New  Discoveries,  Inventions,  or  Patent  Disclosures 

None. 

6  Honors/Awards 

Professional  honors  received  during  career 

A.J.  Laub: 

•  Elected  IEEE  Fellow  in  1986  ("for  contributions  to  algorithms,  numerical  analysis, 
and  mathematical  software  for  control  and  systems  theory") 

•  Distinguished  Member  Award.  IEEE  Control  Systems  Society.  1991 

•  Outstanding  Achievement  in  Teaching  Award.  School  of  Engineering.  University  of 
Southern  California.  May  1983 

•  Academic  Senate  Distinguished  Teaching  Award  (Engineering),  University  of  Califor¬ 
nia.  Santa  Barbara,  June  1993 

•  Control  Systems  Technology  Award.  IEEE  Control  Systems  Society.  1993  ( "for  pio¬ 
neering  efforts  and  continuously  advancing  the  state  of  the  art  in  Computer-Aided 
Control  System  Design" ) 
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